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Motivated by recent density-matrix renormalization group (DMRG) calculations [Yan, Huse, and 
White, Science 332, 1173 (2011)], which claimed that the ground state of the nearest-neighbor spin- 
1/2 Heisenberg antiferromagnet on the kagome lattice geometry is a fully gapped spin liquid with 
numerical signatures of Z2 gauge structure, and a further theoretical work [Lu, Ran, and Lee, Phys. 
Rev. B 83, 224413 (2011)], which gave a classification of all Schwinger-fermion mean-field fully 
symmetric Z2 spin liquids on the kagome lattice, we have thoroughly studied Gutzwiller-projected 
fermionic wave functions by using quantum variational Monte Carlo techniques, hence implementing 
exactly the constraint of one fermion per site. In particular, we investigated the energetics of all 
Z2 candidates (gapped and gapless) that lie in the neighborhood of the energetically competitive 
U(l) gapless spin liquids. By using a state-of-the-art optimization method, we were able to conclu- 
sively show that the U(l) Dirac state is remarkably stable with respect to all Z2 spin liquids in its 
neighborhood, and in particular for opening a gap toward the so-called Z2[0,7r]/3 state, which was 
conjectured to describe the ground state obtained by the DMRG method. Finally, we also consid- 
ered the addition of a small second nearest-neighbor exchange coupling of both antiferromagnetic 
and ferromagnetic type, and obtained similar results, namely, a U(l) Dirac spin-liquid ground state. 



PACS numbers: 75.10.Kt, 75.10.Jm, 75.40.Mg 

Introduction. The nearest-neighbor (NN) spin-1/2 
quantum Heisenberg antiferromagnet (QHAF) on the 
kagome lattice provides ideal conditions for the ampli- 
fication of quantum fluctuations and a consequent sta- 
bilization of an exotic magnetically disordered ground 
state, which may be a valence-bond crystal (VBC)^"^ or 
a spin liquid (SL) with fractionalized excitations.®"^ Re- 
cent experiments have unanimously pointed toward a SL 
behavior;^"^® in particular, Raman spectroscopic data on 
a nearly perfect spin-1/2 kagome compound with Heisen- 
berg couplings (the so-called Herbertsmithite) suggested 
a gapless (algebraic) SL.^^ On the theoretical side, the 
question is still wide open and intensely debated. On 
the one hand, series expansion provided evidence that 
a VBC with a 36-site unit cell has lower energy than 
other proposed competing states.^ On the other hand, it 
was shown that, within the class of Gutzwiller-projected 
fermionic wave functions, a particular algebraic SL, the 
so-called U(l) Dirac state, has a competing energy.^^ Its 
properties were studied in detail in Ref. [19] and it was 
argued that it can be a stable SL state. However, a re- 
cent DMRG study* has challenged the above results, and 
proposed that the ground state can be a fully gapped Z2 
SL with a substantially lower energy as compared to both 
the above estimates. 

The Z2 SLs have the nice property that they are stable 
mean field states and can survive quantum fluctuations. 
Hence, they are more likely to occur as real physical SLs, 
and one can safely use the projective symmetry group 
classification of Z2 SLs beyond mean-field level. This 
complete classification of fully symmetric Z2 SLs on the 
kagome lattice was recently done in Ref. [21] within the 



Schwinger-fermion mean-field theory, resulting in an enu- 
meration of a total of 20 Z2 mean-field states. Their 
main result was the identification of a unique gapped Z2 
SL (called the Z2[0,7r]/3 state) in the neighborhood of 
the U(l) Dirac state. Since the U(l) Dirac SL state has 
the best variational energy among the class of U(l) gap- 
less SLs, in Ref. [21], it has been conjectured that the 
Z2[0,7r]/3 state may describe the ground state that has 
been numerically observed in the DMRG study.* 

In this paper, we thoroughly investigate the possibil- 
ity of any of these Z2 SLs being stabilized as the ground 
state of the NN spin-1/2 QHAF, with a particular em- 
phasis on the Z2[0,7r]/3 state. In practice, we compute 
the energy of optimized variational wave functions that 
are constructed by applying the Gutzwiller projector to 
different states obtained from mean-field Hamiltonians of 
Schwinger fermions. In this respect, by an exact treat- 
ment of the full projector that ensures the one-fermion 
per site constraint, we go much beyond the simple mean- 
field approach of Ref. [21]. We calculate the energies of 
all Z2 SLs which can be realized up to 3rd NN in mean 
field Ansatz and have a non-vanishing 1st NN mean-field 
bond. Only 12 of the 20 Z2 SLs satisfy these criteria, 
and all of them are continuously connected to some U(l) 
gapless SL.^^ Our main result is that, contrary to what 
has been proposed in Ref. [21], the Z2[0,7r]/? state has 
a higher energy than the gapless U(l) Dirac SL, or in 
other words the U(l) Dirac SL is remarkably stable with 
respect to opening of a gap and consequently destabiliz- 
ing into the Z2[0, tt]/? state. We also find that all gapped 
Z2 SLs in the neighborhood of another competing gapless 
state, the uniform resonating-valence bond (RVB) state. 
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FIG. 1. The 1'i^^i^\fi SL ansatz; black (grey) bonds denote 
1st NN real hopping (2nd NN real hopping and real spinon 
pairing) terms; solid (dashed) black bonds have Sij — 1 (—1), 
solid (dashed) grey bonds have Vij — 1 (—1), see Eq. (4). The 
1st NN (2nd NN) mean field Ansatz is written as U^^) — ±0-3 
(U<(,,)> = ±ix2a3 + Aaffi)). The SU(2) flux (P), through 
elementary triangles (e.g., 123) is P123 = as, and that through 
triangles formed by two 1st NN and one 2nd NN bonds (e.g., 
234) is P234 = — (X20"3 + A2cri). Their commutator is non- 
zero, [Pi23,P234] = (— 2ia2)A2. Hence, a finite A2 breaks the 
U(l) gauge structure down to Z2, and opens up an energy 
gap via the Anderson-Higgs mechanism. ^"'^"^ 
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FIG. 2. (Color online) A typical variational Monte Carlo op- 
timization run for the Z2[0, 7r]/3 wave function: (a) variational 
parameters A2, Xz, Mi Cr s^nd (b) energy, as a function 
of SR iterations. In (a), the initialized parameter values are: 
A2 = X2 = 1, A* = -0.8, and Cr = 0.3. The U(l) 2nd NN 
[0,7r;0,7r] Dirac SL corresponds to A2 = 0, X2 = -0.0186(2), 
Cr = 0, as found in Ref. 26. The optimized parameter values 
are obtained by averaging over a much larger number of con- 
verged SR iterations than shown above. In (c), the variation 
in energy upon addition of a small A2 (both for Cr ~ 0, and 
optimized Cr for each value of A2) upon the [0, 7r;0,7r] Dirac 
SL is shown, the increase in energy is apparent. 



have higher energies. Moreover, we find that all Z2 SLs 
have higher energy than the gapless SL states in whose 
neighborhoods they lie. 

Model and wave function. The Hamiltonian for the NN 
spin-1/2 Heisenberg model is 

n = jJ2s.-s,, (1) 

where (ij) denote sums over NN sites and S.; is the spin- 
1/2 operator at site i. All energies will be given in units 
of J. 

The variational wave functions are defined by project- 
ing noncorrelated fermionic states: 

l«'VMc(Xzj,A,,,M,C)) =7'G|*MF(X«j,A,y,/i,C)), (2) 

where Vq = Yli{^~TT'i;t^i,i) ^^^^ Gutzwiller projec- 

tor enforcing the one fermion per site constraint. Here, 
|^MF(Xijj /i, C)) is the ground state of mean-field 
Hamiltonian containing chemical potential, hopping, and 
singlet pairing terms: 

+ J2{{A,,+CS,,)cyi^ + h.c.}, (3) 

where Xij = Xji ^nd Aij — Aji. Besides the chemi- 
cal potential /i, we will also consider real and imaginary 
components of on-site pairing, which are absorbed in 
We briefly mention that a somewhat similar approach, 
based upon a bosonic representation of the spin operators 



(i.e., through Schwinger bosons), has been also used re- 
cently.^^ In the latter case, however, the bosonic nature of 
quasi-particle operators implies that one has to deal with 
permanents instead of determinants, which makes the nu- 
merical calculations much heavier than in our fermionic 
case. 

Different SL phases correspond to different patterns 
of distribution of Xij on the lattice links, along 

with the specification of the on-site terms /i and C- Then, 
a complete specification of a SL state up to nth NN 
amounts to specifying the SU(2) flux through closed 
loops along with the optimized hopping and pairing pa- 
rameters at each geometrical distance. ^"■■^'^ These pa- 
rameters are the Ansdtze of a given state and serve as 
the variational parameters in the physical wave func- 
tion that are optimized within the variational Monte 
Carlo scheme to find the energetically best state. It is 
worth mentioning that we use a sophisticated implemen- 
tation of the stochastic reconfiguration (SR) optimiza- 
tion method, ^^'^^ which allows us to obtain an extremely 
accurate determination of variational parameters. In- 
deed, small energy differences are effectively computed 
by using a correlated sampling, which makes it possible 
to strongly reduce statistical fluctuations. The current 
problem of the study of the instability of a U(l) Dirac 
SL state toward the Z2[0,7r]/3 state will clearly demon- 
strate the power of this method to capture the essential 
subtleties. 

Results. We performed our variational calculations 
on a 432-site cluster with mixed periodic-antiperiodic 
boundary conditions which ensures nondegenerate wave 
functions at half-filling. The large size of the cluster en- 
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sures that the spatial modulations induced in the ob- 
servables by breaking of rotational symmetry (due to 
mixed boundary conditions) remain smaller than the un- 
certainty in the Monte Carlo simulations. 

Among the class of NN fully symmetric and gapless 
SLs, the U(l) Dirac state has the lowest energy. Its en- 
ergy per site is E/J = —0.42863(2), and its Ansatz is 
given by the sign convention for NN bonds in Fig. 1. 
Due to the U(l) flux ip being and tt [exp {iip) = 
ripiaquettc Xij] through triangles and hexagons, respec- 
tively, it is denoted as [0, tt]. Another competing state, 
the NN uniform RVB state has zero flux through any pla- 
quette and is therefore denoted as [0,0]; its energy per 
site is E/J = -0.41216(1). ^^'^^ 

The study in Ref. [21] identified four Z2 SLs in the 
neighborhood of the [0, tt] state; only one of them, the 
Z2[0, 7r]/3 state, was found to be gapped (via the 2nd NN 
spinon pairing term). Its Ansatz up to 2nd NN mean- 
field bond is reproduced in Fig. 1.^^ In a suitable gauge, 
its mean-field Ansatz is specified by five real parameters. 
These parameters are the 1st NN real hopping (xi), 2nd 
NN real hopping (X2), 2nd NN real spinon pairing (A2), 
and two onsite terms, one for the chemical potential /i 
and the other for the real on-site pairing Cr. The mean 
field Hamiltonian can be then conveniently cast in the 
following form: 

■HMF{Z2[0,7r]/3} = Xl ^ Sjj4aCj> 

+ E 4"^^" + ^R(4t4; + '^■c-)}. (4) 

i a 

where {ij) and {{ij)) denote sums over 1st and 2nd NN 
sites, respectively. Sy and Vij encode the sign structure 
of the 1st and 2nd NN bonds, respectively, as shown in 
Fig. 1. The 1st NN real hopping (xi) will be taken as a 
reference, and hence set to unity hereafter. The physical 
variational wave function of this SL state then depends 
on four variational parameters, |^'vMc(X2: Mi Cfl)) — 

^g|*mf(X2, A2,M, Cit))- 

For a generic unbiased starting point in the four- 
dimensional variational space, the variation of parame- 
ters and energy in the SR optimization is shown in Fig. 2. 
As one can clearly see, the energy converges neatly [see 
point B in Fig. 2(b)] to the reference value of the suit- 
ably extended 2nd NN U(l) Dirac SL, the [0, tt; 0, tt] state 
[see point A in Fig. 2(b)], with small but finite X2 [see 
Fig. 2(a)] previously computed by us.^^ For the present 
cluster, these values are E/J = —0.42872(1) per site, 
and X2 = -0.0186(2), = -0.5124(5). Also, it is man- 
ifest that (A2, Cr) ~^ Oj becoming exactly zero (within 
the error bars) after averaging over a sufficient number 
of converged Monte Carlo steps. Here, we bring atten- 
tion to the important fact that, despite the energy having 
converged after sa 400 iterations, the parameters did not 
converge and were still varying, converging to their final 



values much later than the energy (see Fig. 2). This fact 
is possible because, in the energy minimization, forces 
are calculated through the correlated sampling and not 
by energy differences. Our result shows that the en- 
ergy landscape along the manifold connecting the U(l) 
Dirac SL to the Z2[0, 7r]/3 SL is very flat close to the U(l) 
Dirac SL [see Fig. 2(c) for the case Cr{A2}]. Conse- 
quently, a small perturbation around the U(l) Dirac SL, 
e.g., by setting A2 = 0.05 along with the corresponding 
optimized value of Cr = 0.1780(2) will not lead to any de- 
tectable change in energy. Hence, one cannot unambigu- 
ously conclude anything about the stability of the U(l) 
Dirac SL by solely computing the energy of the perturbed 
wave function with fixed parameters, point by point lo- 
cally. Only by performing an accurate SR optimization 
method^'* can one successfully optimize the parameters 
and transparently show that A2 = corresponds to the 
actual minimum of the variational energy. This fact im- 
plies that the U(l) gauge structure is kept intact and the 
Dirac SL state is locally and globally stable with respect 
to destabilizing into the Z2[0,7r]/3 state. We verified this 
result by doing many optimization runs starting from 
different random initial values of the parameters in the 
four-dimensional variational space. 

Regarding the remaining three gapless Z2 SLs in 
the neighborhood of the U(l) Dirac SL, namely, the 
Z2[0,7r]Q;, Z2[0, 7r]7, and Z2[0, 7r](5 states of Ref. 21, our 
study reveals the same result as for the Z2[0, tt]/? state. 
That is, all three of these SLs return back exactly to the 
gapless U(l) Dirac SL state, with the value of the param- 
eter responsible for breaking the U(l) gauge structure 
down to Z2 exactly vanishing upon optimization. Thus, 
we can convincingly conclude that the Z2 neighborhood 
of the U(l) Dirac state does not contain the presumed 
fully gapped Z2 SL found by the DMRG study. 

This conclusion forces us to shift our focus to the Z2 
neighborhood of another fully symmetric (and energeti- 
cally competing) gapless SL, called the uniform RVB or 
the [0, 0] SL. Despite having a slightly higher energy, it 
has the promising feature that all four Z2 SLs in its neigh- 
borhood are gapped, thereby opening up the possibility, 
albeit a slim one, that opening a gap might lead to a large 
gain in energy so as to make one of these four states 
go lower than the U(l) Dirac SL, near to the DMRG 
value of E/J = —0.4379(3) per site. These gapped SLs 
are referred to in Ref. [21] as the Z2[0,0]A, Z2[0,0]B, 
Z2[0,0]C, and Z2[0,0]D states; for their Ansdtze, see Ta- 
ble I in Ref. [21] and also the supplementary material. 
Our simulations show that all four of these SLs return 
upon optimization to the gapless uniform RVB SL, with 
optimized Xn- In particular, case by case we see that: for 
the Z2[0, 0]A SL, the 2nd NN spinon pairing term goes to 
zero along with the on-site pairing term, thus returning 
back to the 2nd NN uniform RVB SL, the [0, 0; tt, tt] state 
given by optimized X2 = -0.032(1);26 for the Z2[0,0]B 
SL, the NN spinon pairing term goes to zero upon op- 
timization, thus giving back the NN uniform RVB SL; 
the Z2[0, 0]C SL upon optimizing flows to the 3rd NN 
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uniform RVB SL with optimized xi ^-^id X3S, with the 
spinon pairing term at 3rd NN becoming zero; and fi- 
nally, the Z2[0,0]D SL flows back to the 2nd NN uni- 
form RVB SL. The results showing how the energies of 
these extended gapless uniform RVB SLs increase as the 
U(l) —7- Z2 gauge breaking parameter is tuned on from 
zero to a small finite value are reported in the supple- 
mentary material (at the end of the article) . 

For reasons of completeness, we mention that there 
are two more gapless U(l) SLs in whose neighborhoods 
the remaining four %i SLs (all gapless) lie.^^ However, 
these U(l) SLs suffer from a macroscopic degeneracy at 
half-filling which leads to an open shell. This degeneracy 
being macroscopic cannot be removed by using any of the 
four real boundary condition possibilities. Hence, their 
energy can only be computed approximately in the limit 
by inserting an additional 6 fiux through the triangle mo- 
tifs and consequently removing n ~ 29 through hexagon 
motifs, and then taking the limit 6 —> 0. The energy of 
the SL-[7r, tt] computed in this way is E/J ~ -0.38372(1) 
per site, which is much higher than those of other gapless 
U(l) states. Hence, we did not carry out an analysis of 
Z2 SLs in these two neighborhoods. 

Finally, we also investigated the possibility of stabiliza- 
tion of the Z2[0,7r]^ state upon addition of a small 2nd 
NN exchange coupling (J') of both antiferromagnetic and 
ferromagnetic type in the NN spin- 1/2 QHAF. In both 
cases, on optimization we found that (A2, Cr) ~^ 0, be- 
coming exactly zero (within error bars) after averaging 
over a sufficient number of converged Monte Carlo iter- 
ations. Thus, we recover the suitably extended, gapless 
2nd NN U(l) Dirac SL. In particular, for J' /J = 0.10, 
this is the [0, tt; tt, 0] state with optimized X2 — 0.0924(2), 



and E/J = -0.43200(2) per site; for J' /J = -0.10, this 
is the [0, 7r;0,7r] state with optimized X2 — —0.1066(2), 
and E/J = -0.42898(2).26 

In summary, we investigated the possibility of stabiliz- 
ing gapped Z2 SLs in the NN and next-nearest-neighbor 
(NNN) spin-1/2 QHAF on a kagome lattice. We found 
that none of the five gapped Z2 SLs [one connected to 
the U(l) Dirac state and the other four connected to 
the uniform RVB state] can occur as ground states. In 
particular, the most promising gapped SL conjectured to 
describe the ground state, the Z2[0,7r]/3 state, is always 
higher in energy than the U(l) Dirac SL. Our system- 
atic numerical results bring us to the conclusion that, at 
least within the Schwinger fermion approach of the spin 
model, the U(l) Dirac SL has the best variational en- 
ergy for the NN and NNN spin-1/2 QHAF on kagome 
lattice. The confiict of our results, which point toward a 
gapless ground state, and the ones by recent DMRG cal- 
culations, which instead suggested the presence of a fully 
gapped spectrum, remains open and deserves further in- 
vestigations. One possible direction would be to consider 
further improvements of our variational wave functions, 
based upon the application of few Lanczos steps or an 
approximated (fixed-node) projection technique, which, 
e.g., gives an energy of E/J — —0.431453(2) per site for 
the NN U(l) Dirac SL, and E/J = -0.431443(2) for the 
NNN [0, tt; 0, tt] state. Another possible direction would 
be to explore the energetics of gapped Z2 SLs which break 
some symmetries such as time-reversal. The possibility 
that the fully gapped SL found by the DMRG study pos- 
sesses a different low energy gauge structure other than 
Z2 also remains open. 
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In this supplementary material, we present the explicit plots showing the variation of energy 
as a function of the (7(1) — > Z2 gauge breaking parameters for the four gapped Z2 SLs in the 
neighborhood of the gapless [0,0] SL (the uniform RVB state). We also reproduce the ansatz (from 
Table I of Ref. 21) of the five SLs investigated in paper, so as to make the paper self contained. 
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TABLE I. The mean field ansatz of tlie five gapped SLs investigated by us, given only up to the neighbor at which the gauge 
symmetry is broken, in a form used by us in numerical simulations. The parameters highlighted in red are responsible for 
opening a gap by breaking the (7(1) gauge symmetry down to Z2. The Usrdn.n. denotes bonds of length 2 connecting two sites 
and passing through a third site (such as the bond 1 — >■ 4 in Fig. 1 in the main paper); instead, U^rdn.n. denotes bonds of length 
2 which dont pass through any site. 
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FIG. 1. For the four Z2 spin liquids in the neighborhood of the [0,0] state: The manner of variation of energy as the gauge 
breaking parameter (shown in red in Table I) is tuned on from zero to a small finite value, is shown. The parameters in black 



in Table I are fixed to their optimized values, and correspond to a suitably extended n 
increase in energy upon opening a gap is apparent. 



NN gapless uniform RVB SL. The 



In the main paper, we have discussed the numerical results of a full Monte Carlo optimization for the four gapped 
Z2 SLs in the neighborhood of the [0, 0] state. In particular, we find that, by starting from arbitrary points in 
the variational space, one returns back to the gapless [0,0] reference SL (suitably extended to n^^ NN), with the 
value of the gauge breaking parameter going to zero. The energies of the NN [0, 0] SL and the extended 2nd and 
3rd n.n. [0,0] SLs correspond to the points B, A, C in Fig. 1 respectively, these energies are E/J = —0.41216(1), 
E/J — —0.41228(1), E/J — —0.412308(1) respectively. Here, we show the local explicit manner of increase in energy 
of these extended gapless [0, 0] SLs upon addition of a small gap opening (gauge breaking) parameter (given in red in 
Table. I), like we did for the Z2[0, 7r]/3 state in Fig. 2(c) of main paper. It should be noted that the case of Z2[0, 0]A 
with Cr = is equivalent to the Z2[0,0]D case. Hence, we have the same plot for both SLs in Fig. 1. 



